Well-tempered metadynamics: a smoothly-converging and tunable free-energy method 
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We present a method for determining the free energy dependence on a selected number of col- 
lective variables using an adaptive bias. The formalism provides a unified description which has 
metadynamics and canonical sampling as limiting cases. Convergence and errors can be rigorously 
and easily controlled. The parameters of the simulation can be tuned so as to focus the computa- 
tional effort only on the physically relevant regions of the order parameter space. The algorithm is 
tested on the reconstruction of alanine dipeptide free energy landscape. 
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Computing free energy differences is of crucial impor- 
tance in molecular dynamics (MD) and Monte Carlo 
(MC) simulations. Whenever it is possible to define a few 
collective variables (CVs) that provide a coarse-grained 
description of the slow modes [2, HJ, it is also of great 
relevance to compute the associated free energy surface 
(FES). In order to draw such a surface a straightfor- 
ward approach is often not possible due to high barriers 
or other sampling bottlenecks. A standard strategy for 
overcoming this problem is to introduce an external bi- 
asing potential that forces the system to explore regions 
of high free energy Q. A major progress has been the 
recent introduction of adaptive non-equilibrium methods 
0, H, H, 0] ■ In all these methods the simulation history is 
used to enhance the sampling speed. In a MC run. this 
can be done by varying the MC acceptance probability 
every time a new configuration is visited 0] , while in MD 
a time-dependent bias can be added either to the force 
[|| or to the potential @, [| . 

We shall focus here on metadynamics @, which has 
proven its effectiveness in a variety of contexts [E [E EE 
El EE EE EI EE EE EE El ■ In metadynamics the sys- 
tem evolution is biased by a history-dependent potential 
that is constructed as the sum of Gaussian functions [l9[ 
deposited along the trajectory in the CVs space. After a 
transient, the bias potential compensates the underlying 
FES and provides an estimate of its dependence on the 
CVs. A formal justification of this procedure has been 
given in Ref . [IE] ■ In spite of its success there is a need to 
improve metadynamics in several respects. First of all, it 
is often difficult to decide when to terminate a metady- 
namics run. In fact, in a single run, the free energy does 
not converge to a definite value but fluctuates around the 
correct result, leading to an average error which is pro- 
portional to the square root of the bias potential deposi- 
tion rate [IE HH . Reducing this rate implies increasing 
the time required to fill the FES. Furthermore, in prac- 
tical application, continuing a run carries the risk that 
the system is irreversibly pushed in regions of configu- 



rational space which are not physically relevant. These 
issues have already been recognized and different ad-hoc 
solutions have been proposed to alleviate these problems 

[EI EE HE HE HI- 

In this Letter, inspired by the self-healing umbrella 
sampling method 0, we substantially improve metady- 
namics such that we obtain an estimate of the FES that 
converges to the exact result in the long time limit. Con- 
trary to ordinary metadynamics, our approach offers the 
possibility of controlling the regions of FES that are phys- 
ically meaningful to explore. Besides being highly effec- 
tive and controllable this new method provides a uni- 
fied framework whose limiting cases are standard meta- 
dynamics and non-biased standard sampling. We dub 
this new scheme well-tempered metadynamics. 

Let us consider a system described by a set of micro- 
scopic coordinates q and a potential energy U(q), evolv- 
ing under the action of a dynamics (e.g. MD or MC) 
whose equilibrium distribution is canonical at the tem- 
perature T. We want to determine the free energy de- 
pendence on a set of collective variable s(q). The FES 
can be written within an immaterial constant as 



F(s) 



-T lim In iV(s,t), 

t — >-00 



(1) 



where N(s, t) — j Q S s s ^ t i^dt' is the histogram of the vari- 
able s obtained from an unbiased simulation. By con- 
struction, N(s,0) = and its time derivative N(s,t) = 



»(*)• 



To accelerate sampling we bias the dynamics by 



adding the history-dependent potential 



wN(s,t) 
AT 



V(s,t) = AT In 1 



(2) 
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where w has the dimension of an energy rate, AT is a 
temperature and N(s,t) comes from the biased simula- 
tion. Since V is a monotonic function of A, such a bias 
potential disfavors the more frequently visited configura- 
tions. A crucial quantity is the rate at which the potential 
is modified. In particular, slower variation rates lead to 
a dynamics of the microscopic variables q which is closer 
to thermodynamic equilibrium. From Eq. ^ it follows 
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that the rate with which V(s, t) changes is: 



V(s,t) 



ujATS 



s,s(t) 



AT + LuN(s,t) 



cue 
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(3) 



The connection with metadynamics is evident if we ex- 
amine Eq. ([3]) and replace S s ,s(t) with a finite width Gaus- 
sian. Therefore, our scheme can easily be implemented 
in any metadynamics code by rescaling the height of the 
Gaussians according to Eq. ([3]). Using the notation in 
Ref. [2l[, the height of each Gaussian is determined by 

_ V*(s,t) 

w = toe at TG; where tq is the time interval at which 
Gaussians arc deposited. Thus to represents the initial 
bias deposition rate. 

Two important properties need to be underlined. The 
first is that since the histogram N(s, t) grows linearly 
with simulation time, the rate V(s,t) tends to zero as 
oc 1/t. This is the simplest, if possibly not the optimal 
[2f| , way to have a rate decrease fast enough for the bias 
eventually to converge, yet slow enough for the final re- 
sult not to depend on the initial condition V(s,0). Sim- 
ilar arguments have been used in the field of stochastic 
optimization [26|, [27j . The second property is that V is 
not uniform in the s space since at a given point the rate 
is inversely proportional to the time already spent there. 
This latter feature distinguishes our approach from oth- 
ers in which 1/t strategy has also been suggested either 
explicitly [H| or implicitly Q. 

For large times, V(s, t) varies so slowly that one can 
assume that the q's reach equilibrium, the probability 
distribution becomes P{s, t)ds oc exp y F ( s )^ v ( s ^) ^ s 
and one has: 



V(s,t) = we zr^P(s,t) = toe ar" 



F(s) + V(s.t) 



Ids 

(4) 

This implies that V(s,t — > oo) ~ — AT+T F(s), modulo 
a constant. Thus at variance with metadynamics and 
other methods the bias does not fully compensate F(s), 
rather one has that F(s) + V(s) = T ^ AT F(s) leading to 
the following distribution of s: 

P(s,t — > oo)ds cx e T+ATds. (5) 
In practice, using Eq. ([2]) the FES can be estimated as 



F(s,t) 



T + AT 
AT 



V(s,t) = -(T+AT) In 1 



u>N(s,t) 



AT 

(6) 

Let us examine the two limiting cases, AT = and 
AT — > oo . For AT ~ the bias is equal to zero 
and Eq. ([6|) reduces to Eq. ([1]). More interesting is the 
AT — > oo limit. In this case, the deposition rate is con- 
stant, and from Eq. © one finds that F(s, t) = — V(s, t) 
and the standard metadynamics algorithm is recovered. 
Note however that the limit AT — > oo is singular: if we 
first let AT — > oo, the convergence of V(s,t) for t — > oo 
cannot be demonstrated by means of Eq. This is a 



reflection of the already noted drawback of metadynam- 
ics that in a single simulation, the bias does not converge 
but oscillates around the correct F(s) value. In interme- 
diate cases the calculated FES is the one corresponding 
to the target temperature T, with the transverse degrees 
of freedom correctly sampled. However, the s probabil- 
ity distribution is altered and corresponds to an enhanced 
temperature T + AT. It must be stressed that this re- 
sult has been obtained without having to assume adia- 
batic separation between s and the other variables as in 

Rcfs. [la, [S3, EH. 

Much is to be gained computationally by well- 
tempered metadynamics. By tuning AT one can increase 
barrier crossing and facilitate the exploration in the CVs 
space. Furthermore using a finite value of AT one auto- 
matically limits the exploration of the FES region to an 
energy range of the order T + AT. Hence the exploration 
of the FES can be limited to the physically interesting 
regions of s. Longer simulation time results in improved 
statistical accuracy in the relevant regions. The risk of 
overfilling is avoided and optimal use is made of the com- 
puter time. Deciding when to stop the run is now simple 
and post-processing [1, [22| is not necessary. 

As an illustration we study the FES of alanine dipep- 
tide in vacuum as a function of the backbone dihedral 
angles ($,^). This surface has been well studied and is 
known to exhibit two minima CV eo and Cj ax separated 
by a barrier of ~9 kcal mol -1 [32|, |33|. Since such a bar- 
rier cannot be crossed with standard dynamics at room 
temperature this system has provided a testing ground 
for many sampling schemes. The CHARMM27 [34| force 
field has been used in ORAC MD code [35[ and canon- 
ical sampling at a temperature of 300 K was achieved 
by means of the stochastic thermostat in Ref. pBtij . The 
Gaussian width was set to 20 degrees, and the deposi- 
tion interval was 120 fs with a starting Gaussian height 
of 0.287 kcal mol -1 , which corresponds to a deposition 
rate w=2.4 cal mol _1 fs -1 . 

We calculated a reference "J) using standard um- 
brella sampling which is in good agreement with previ- 
ous studies. On this surface we superimpose three dif- 
ferent trajectories (see Fig. [T]) started from the same ini- 
tial conditions \pieq (-83,74)], but with three different 
choices of AT (600 K, 1800 K and 4200 K). In all three 
cases the secondary metastable state C-j ax — (70, —70) 
was frequently visited. It is worth noting that, as dis- 
cussed earlier, by increasing AT larger and larger re- 
gions were explored. In order to demonstrate how the 
method converges, for the three mentioned cases, in 
Fig. [1] we also show the time evolution of AF(t) — 
F(C7 aX7 t) — F(Cj eqi t), i. e. the estimated free energy 
difference between the two minima. AF(t) converges to 
the reference value (AT «2.2 kcal mol -1 ) in all three 
trajectories. At variance with standard metadynamics, 
the time derivative of the bias potential tends to zero 
and the fluctuations around the correct value are pro- 
gressively damped. All three simulations provide an ac- 
curate estimate of the free energy difference within a few 
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FIG. 1: Panels (a, b, c): Green dots represent 6 ns long 
trajectories in the ($, space for different choices of AT 
[600 K(a), 1800 K (b), and 4200 K (c)]. The underlying color 
map (kcal mol -1 ) shows the reference free energy landscape. 
Panel (d): Estimate of the free energy difference between the 
two metastable minima C-t a x (70,-70) and CVe, (-83,74) as a 
function of the simulation time, as obtained from the same 
trajectories. 
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FIG. 2: Panel (a): Time evolution of (e(t))V* for differ- 
ent choices of AT. (e(t)} is the error as defined in Eq. 0, 
averaged over an ensemble of 100 independent atomistic sim- 
ulations starting from Creq- Panel (b): Dependence of k (see 
text for definition) on AT, as estimated from 6 ns long tra- 
jectories. 



nanoseconds, even in the lowest AT case where the lower 
number of barrier crossing events leads to a jumpier AT 
evolution. 

As a measure of the error of F(^, <£) in the relevant 
regions, we define 

e(f) ={^J $ ) - ^(*. $ ' *) - C{t)fd<S>d^J ' 

(7) 

where T is the region in dihedral space such that 
(T(*,$) - F(C 7eq )) < 10 kcal, and A is its area. V 
is defined to include all the minima and all the transi- 
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FIG. 3: (color online) Time evolution of the average error 
(e(i)} for different values of tb and D*, where D$ = 12.3 
deg 2 fs _1 and AT=1200 K. The error is averaged over an 
ensemble of 1000 independent Langevin simulations starting 

from Cjeq- 



tion states. The value of C(t) is chosen so as to align 
the averages of F and T over T. It is seen that after an 
initial transient period, (e(t)) converges to zero as k/y/i. 
Such behavior is shown in Fig. HJa) where {e{t))\/t is 
plotted against the simulation time for three values of 
AT. This is clearly at variance with standard metady- 
namics in which the error does not to converge to zero 
during a single simulation 0, [2l[ . The behavior of the 
present scheme is consistent with an error analysis done 
on a simulation performed at a constant bias. In Fig.[5{b) 
we study the dependence of k = limt_ >00 (e(f))vt on AT 
as a way of optimizing the choice of AT. In this case 
the optimal choice is close to AT = 1200/r resulting 
in a sampling temperature for the collective variables of 
T + AT = 1500^, which is of the order of magnitude of 
the barrier height. Its actual value may depend on the s 
relaxation times and on the area one wishes to explore. 

We discuss now the role of ui, the initial deposition rate 
which we relate to the time constant tr = that sets 
the time scale for the bias evolution. While in the long 
time limit tb is irrelevant, it could affect the transient 
regime in a non-trivial way. At constant AT, a small 
tb implies a high initial deposition rate, thus leading to 
rapid filling of the wells. However, if tb is too small 
relative to the time necessary to properly average out 
the transverse degrees of freedom, the large fluctuations 
in the initial FES reconstruction need a longer time to 
be recovered. 

This effect is conveniently investigated by introducing 
an artificial model based on the alanine dipeptide FES. 
We model the dynamics on the two-dimensional space 
((f>, 4') with a high-friction Langevin equation driven by 
the free energy surface F ($, >F) and the diffusion coef- 
ficients D^, determined from the atomistic simula- 
tions. We shall apply our scheme to calculate the one- 
dimensional projection F(<&), using a one-dimensional 
bias on In such a case the relaxation speed of the 
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transverse degree of freedom ^ can be tuned by changing 
D\j,, thus mimicking a situation in which the transverse 
degrees of freedom are either fast or slow. As can be 
seen in Fig. [3J in the fast case, the orthogonal degree of 
freedom is rapidly averaged out, resulting in a Markovian 
dynamics on ^, and a small tb is the best choice. In the 
slow case, the effective dynamics of <3> is strongly non- 
Mar kovian due to coupling with 'J, and a small tb is not 
the best choice since it results in an increase of the tran- 
sient time. However, it is worth noting that the method 
is robust and in the range of reported cases, which spans 
two orders of magnitude in tb and D^ !l the calculation 
converged to the same results on approximately the same 
time scale. 



In conclusion, well-tempered mctadynamics solves the 
convergence problems of mctadynamics and allows the 
computational effort to be focused on the physically rel- 
evant regions of the conformational space. The latter 
property makes it possible to use adaptive-bias methods 
in higher dimensionality cases, thus paving the way for 
the study of complex systems where it is difficult to se- 
lect a priori a very small number of relevant degrees of 
freedom. The proposed approach can easily be applied to 
generalizations of mctadynamics based on multiple repli- 
cas [l5l.[37lj38j. and can be extended to the Wang-Landau 
algorithm [4|. 
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